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CO ' Periodic boundary conditions are widely used in the simulation of systems with an extremely 

\ large number of particles, and the period vectors are the degrees of freedom replacing those of the 

. image particles. We derived dynamical equations for the periods by applying Newton's Second Law 

to halves of the macroscopic system with external forces considered explicitly in the case of constant 
I stress. The resulting internal stress has both the interaction term and the controversial kinetic- 

^ , energy term. A new, statistical explanation for the interaction term was given. The kinetic-energy 

' term was obtained by considering collisions between particles and walls, as well as introducing an 

additional "force" associated with the pure transport of momentum. For assumed fixed periods at 
04 , low temperature, it can be used to calculate the external stress under which the system is in an 

equilibrium state. We gave an example calculation for Cobalt crystal structures under constant 
stress. 
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I. INTRODUCTION 



With Newtonian dynamics, computer simulations of particle systems play an important role in Condensed Matter 
Physics, Biophysics, Chemical Physics, Space Physics, and Astrophysics. Due to the extremely large number of 
particles, periodic boundary conditions (PBCs) are often employed. So the whole system is treated as a crystal 
composed of periodic cells, i.e., from the condensed matter physics point of view. Since the PBCs are imaginary 
conditions, are we allowed to choose the periods or the cell edge vectors arbitrarily, even for free systems at low 
temperature? We know that there is a great variety of period values in real crystals. So PBCs are acceptable in some 
Q ■ sense. But crystal periods are not arbitrary. Neither should be the periods in PBCs. The crystal periods are also 
' called lattice constants, which are not changing in low temperature equilibrium states, under given constant external 
forces or pressure (normal force per unit area) or stress (force per unit area), which act on the surfaces of crystals 
only. However if these external forces are subject to change, the periods will definitely change accordingly. In fact, the 
periods were already recognized as new degrees of freedom replacing those of the image particles previously P|. This 
fS| is not difficult to understand. If the positions of all particles in one cell are determined, those of particles in other 
C — ■ cells are also determined if and only if the periods are determined. So we need to develope dynamical equations for 
I the periods with consideration of external stress, rather than assigning some estimated fixed values to them. Another 
related problem is how external forces affect internal particles. If the periods are affected by the external forces acting 
on particles at the boundaries, as most of the internal particles' positions are also determined by the periods, the 
external forces eventually affect the internalparticles' behavior. 

As molecular dynamics (MD) simulation 0, 0, 01 can be used to simulate all of the above processes in various fields 
of physics, let us use the MD terminology in the following. The center cell is called "MD cell" , and the particles in 
g ; the MD cell are called "MD particles" . 

I • The first periodic dynamical equation was the one for cell volume proposed by Andersen in 1980j5(| in his simulation 
' O \ of fluids subjected to a constant external pressure. Shortly thereafter, Parrinello and Rahman extended Andersen's 
1 idea to a new MD theory (PRMD), in which the dynamical equations for the complete periods were proposed for 
constant external pressure and constant general external stress respectively pj, i6|. Since then, many phase transition 
simulations on crystal structures induced by external forces have been carried out[3, H, 0, llJj especially after it 
was combined with the well-known Car-Parrinello MDfl^ ll^ E^. Later, some various versions of the PRMD were 
also proposed 15, 16, 17, .18, ,19, .2(I.2h|. 

However one problem in these works is the violation of Newton's Second Lawj:22|j. They all use the scaled particle 
' position vectors = (^i,?7i,Ci) defined in = £^iSi+rjib+C,iC = Hs^, where is the MD particle i's position vector, a, 
b and c are the cell edge (period) vectors (forming a right-handed triad), and H is the matrix formed by {a, b, c}, 

as generalized coordinates. By taking (Hs^)^, rather than ^Hsj -f- Hs^^ , as kinetic energy of particle i 

with mass in the Lagrangian of the MD cell and applying the Lagrangian dynamics, they obtained dynamical 
equations for both the MD particles and the periods. Although these deduction can lead to the generally believed 
kinetic-energy term of the internal stress and interaction term of the internal stress, unfortunately, the generated 
dynamical equations for the MD particles from the incomplete kinetic energy violate Newton's Second Law[22. If the 
complete particle kinetic energy were used in the Lagrangian of their works, the kinetic- energy term of the internal 
stress would have not been generated in physics, as in our previous works |2a 12^. 
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As a matter of fact, using s^, rather than r^, as generaUzed coordinates, can not simphfy but comphcate the formula, 
and is unnecessary to couple the particle positions and the periods , ,23i. .24] . The particle positions and the periods are 
coupled when particles of more than one cell are considered. This paper adopts the MD particle positions and the 
periods, the natural and complete degrees of freedom, as independent dynamical variables. 

As external forces act on particles in surfaces directly only, it is usually difficult to present a clear physics picture 
for their action in a Lagrangian of the MD cell. So in this paper, we will apply Newton's Second Law to halves of the 
macroscopic crystals to get dynamical equations for periods, with explicit action of the external forces on the surfaces 
of the crystals, and keep Newton's Second Law for the MD particles in the original form to avoid the above violation 
problem. By also considering the external forces involved in collisions between infinitely far away walls and particles 
near the walls and introducing an additional "force" accompanying pure transport of momentum, the controversial [2^ 
kinetic-energy term of the internal stress is found in this paper. 

This paper is organized as follows. The model and basic ideas are given in Sec. ^ The simplest case of low 
temperature without collisions is discussed in Sec. IIIII then a new statistical explanation of the interaction term is 
reported in Sec. IIVI The collision between particles and walls is discussed in Sec. |3 which results in the kinetic- 
energy term of the internal stress eventually. A further discussion on the kinetic-energy term is presented, in which 
the additional "force" in pure transport of momentum is introduced, in Sec. IVII Calculated cobalt structures under 
constant stress is shown in Sec. I VIII and Sec. IVIIII is devoted to a summary. 



II. MODEL AND BASIC IDEAS 



The macroscopic bulk of a material with a periodic structure expressed by H is taken as the model for the study. 
For simplicity, let us suppose that it consists of complete cells only, denoted by the corresponding lattice translation 
vectors T = T^a. + T^h + TcC, with Ta, Tb, and Tc integers, and T = for the MD cell. We study the properties 
of the inner part of the bulk around the center, from which the surface can be regarded infinitely far away, and only 
averaged surface effects affect the inner part. 

The external action on its surface is expressed by the constant external stress tensor (or dyad) T • It can be 
divided into two parts. One is the collision interaction between particles and walls surrounding the bulk, denoted by 

the stress tensor C , and the other is the rest external interaction, denoted by the stress tensor P . The former is 
particle-velocity related and usually not significant at low temperatures. The following section will address the case 
with the latter considered only. Later, the whole external interaction 

T = C + P (1) 

will be dealt with. 

By definition, the external force acting on an infinitesimal surface area vector ds of the bulk is 

dF = ^ -ds. (2) 
Here, the direction of ds is from inside to the outside of the bulk. The net external force on the bulk is 

T -ds = T ■ f ds = 0, (3) 

surface J surface 

i.e., there is no acceleration for the mass center of the bulk. Let us set the center of mass at a fixed position. Usually, 

the external stress T is assumed to be symmetric, i.e., for all of its components Ti j = Tj i. This stress symmetry 
ensures that the net external torque on the bulk is zero. This can easily be seen in a case of a cuboid, whose six 
surfaces are in the directions of the unit axes ie^, i^y) ^^'^ of a Cartesian coordinate system with a volume 

V. Then the total torque in the direction is (T^:,;^ — Tj,^^;) V for T acting on the six surfaces of the cuboid. We 

further assume that the bulk does not rotate either. However, under the constant external stress T , the macroscopic 
volume and shape of the bulk may change, or the bulk may expand or contract with time, as it may not be in a static 
equilibrium state. Of course, the microscopic positions of the MD particles and the volume = (a x b) -c and shape 
of the MD cell also change correspondingly. 

The dynamics for the MD particles is well known. Let us repeat Newton's Second Law for these particles: 



rriiri = F^, (i = 1,2,---, n) , 



(4) 
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where n is the number of particles in each cell and is the net force acting on MD particle i from all other particles 
without external force, which does not act on the MD particles directly. The corresponding equilibrium condition for 
these degrees of freedom is 

F,; = 0, (^ = 1,2,..., n) . (5) 

Simplifying further, we only discuss pair potentials ip{ri j) — (p{\ri — rj|). The extension to many-body interactions 
is straightforward by examining our previous work 24J. It is easy to deduce that 

n n 

^F, = = ^m,r„ (6) 

i=l 4=1 

with the help of Newton's Third Law for interactions between particles, and the periodicity of the crystal. 

For the bulk, the periods H arc the only left degrees of freedom. In order to find dynamical equations for them, we 
use a plane PP' to cut the bulk into a right part and a left part, with S as the area vector of the cross section between 
the two parts in the direction of pointing the right part, as shown in Fig. 1. Plane PP' is parallel to CTh = dfl/dh, 
the surface area vector of a cell, with h = a, b, or c. The right (i?) part contains all and only T cells with Th > 0, 
and the left (L) part contains all the rest cells. Applying Newton's Second Law on the R half crystal is the starting 
point in seeking the dynamic equations for periods, while keeping Eq. Q in the original form. 



III. LOW TEMPERATURE DYNAMICS 



For low temperature, we can ignore collisions between particles and walls, which will be addressed in the next 
section. Let Fl^h be the net force of the L part acting on the R part. By using the net external force acting on the 
R part 



Fe,r = P ■ds = P ■ ds= P -S, (7) 

jR,sf J R,sf 

where ^.i is over the surface of the bulk in the R part, the dynamical equation of the R part is 



MrTrc =Yl^r+ P -S, (8) 



where Mr and yrc are the total mass and acceleration of the center of mass of the R part. The corresponding 
macroscopic equilibrium condition for the R part is 



' L^R- 



P -8 = 0. (9) 



The R part can be in an equilibrium state, with Eqs. Q and jsj satisfied and the net torque zero, if Y^^r is 
uniformly distributed across the section S. Let us assume this uniformity for Y^^r and divide Eq. © by 

A^h = ISI/khl , (10) 

then 

1 

— MrTrc = Fh+ P -CTh , (11) 

where Fh is the net force of the L part acting on the half-line-cell bar Bh composed of the MD cell, cells h, 2h, 3h, 
. . ., up to the right surface, as shown in Fig. 1. It is easy to obtain 

(Th<0) „ (Th>0) „ ^ n 

Fh = X! X! ^hfi,0^i,T = X! X! ?h£i,0^j,T = 2 X! X! ^hfi,0-+j,T, (12) 
T = l T = l T^ai,j = l 

where fi.x^j.T' is the force acting on particle j of the T' cell by particle i in the T cell|23|. 
With the cell potential energy of interactions among particles 

n n 

E= ^ ^(|r,-r,|) + -^^^(|r,-r,-T|), (13) 

i>j = l i,i = lT^O 
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let us define a main interaction tensor (dyad) as 



dE\ f dE 
a+ 



9a 



dh 



b+ 



dE 
dc 



which can be easily educed to the form 



We then have 



2n 



Fh— £m -fh- 



(14) 



(15) 



(16) 



Using Eq. © and the assumption that there are complete cells only in the bulk, the left hand side of Eq. Hll|) 
becomes 



TS-R i=l TeR 

where Mceii = ™i- Noticing T = Taa + Tbb + Tc c, and introducing 



Th 



nNh 



TeR 



(17) 

(18) 

(19) 
(20) 
(21) 

(22) 



with the matrixes T = {Ta, Tb, Tc}, cr = {(Ta, <^h, o'c}- The corresponding equilibrium condition for the H 
degrees of freedom becomes 

em + ^P^ = 0, (23) 
= fl^ ^ 0. This equilibrium condition is equivalent to the macroscopic equilibrium condition of Eq. ll^. 



as well as defining the dyad 

D = aCTa + bCTb + CCTc, 

we obtain 

Gh ^'d^ Th. 

Then Eq. ifTTl) becomes 

D -Th = (^Ern + P ^ • CTh, (h = a, b, or c), 
which may be written as the dot product 

D ■ T = f £m + p)-^o^ , 



as 



Let us further expand the Eq. H17|l to 



where 



Gh = ah,aa + ah,bb + Q!h,cC, 



(24) 



(25) 



Te-R 



As the summation of the positive values of Th'^^th is nearly the same as the absolute summation of the negative values 
of Th'^h, and Th is always positive in the R part, we can make the assumption ah,h 3> |ah.h'7th|- Then Eq. H22|) 
becomes 



H' = e™ + P ]■ a , 



(26) 
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with the matrix H' — |Q:a.aa-, Q^b.bb, ac.ccj-. Considering that ah,h is determined by the size and shape of the 
macroscopic bulk, and that we are interested in properties of the inner part of the bulk, which should be unaffected 
by them, we can further set aa,a = ctb.h = ctc.c = W . This may be regarded as a support for the fictitious kinetic 
energy for the periods in the PRMD. Then Eq. (|26|l further becomes 

WU = ( £^ + P^) ■ . (27) 

For the case of constant external pressure p , 

T= -p T , (28) 

with identity tensor / , then Eq. (|27|l tm'ns into 

WH = (^tn -p) ■ ^<y^ , (29) 
which is the same as the dynamical equation in our previous work [2^. 



IV. INTERACTION TENSOR 

Eq. suggests that the main interaction tensor Em can be regarded as internal stress, which is usually believed 
to have a kinetic-energy term and an interaction term. Zhou proved that the kinetic-energy term has nothing to do 
with mechanical or regular force or stress very recentlyp^. However we will find it by considering collisions between 
walls and particles, and by introducing a non-regular force later. Here let us focus on the interaction term. It usually 
takes the form of 

i>j 

with particles / = (i, T) and J — (j, T') running over all the material and V being the whole volume of the material. 
Haile provided a derivation of Eq. H30|l by taking an unweighted average of internal pressure or stress components 
over all possible locations of parallel cross sections of a material in Appendix B.3 of his bookQ. Using Fig. 1, let 
us do a similar unweighted average but of the dynamical equation, Eq. ((HJ, over all possible locations of a plane 
QQ' passing through the MD cell and parallel to the plane PP'\\cr-^, supposing the cross section area vector S of the 
crystal in QQ' remains unchanged. The particles' availability and configuration surrounding a surface particle and 
those surrounding an inner particle are different. For simplicity, we use those surrounding inner particles for surface 
particles, as the contribution to Eq. H3U|) from the surface is apparently much smaller than that from the inner part. 

For a given location of the plane QQ' , all particles on its left should be excluded from the left side of Eq. ||SJ). The 
total probability for the MD particle i to be excluded from the left-hand side of Eq. ((SJ is (h — r^) ■ a\i/Vl, also the 
total probability for the plane QQ' appearing on the right-hand side of the particle. So the averaged Eq. © over all 
QQ's is 

mrvrc -n^Y. rf-^"'^'''^ (3^) 

i 

By considering the periodicity of the bulk, Eq. (|3U|) can be simplified as 

/ n n \ 

= ^ E (f.-^o-'.o) (r. - r,) + E E (f^-^T-,o) (r. - r, - T) 

1 / " " \ 

^ 2Vt\^ (fi,o^»,or» fj,o^i,orj) + E E (f?>T-*,ori + fi,-T^j,orj -I- f,,o^j,TT) 

T#0»,i = l / 
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1 " 



T#0 i,j = l 



1 

m 



with 



1 " 



r2 



(32) 



(33) 



1=1 



With Eqs. igj, ©, ((Tn|l . lfTC|) and used, Eq. ||2U reduces to Eq. ©. Consequently, the dynamical equation 
remains unchanged even when the average over all QQ's is carried out, and the main interaction tensor contains 
the main part of the interaction term. It is also deduced that if the bulk is in an equilibrium state with both Eqs. Q 
and H23|l satisfied for a given plane PP' , Eq. Q is also satisfied with respect to any location of QQ' . 

One common mistake is the misuse of as £ , with £,„ unrealized. In fact , Sm is the major part of the 
interaction term. For periodic systems, Eq. H14() can be easily extended to the cases for many-body potential or 
potential calculated from the quantum mechanics, while Eq. H3U|I may not be so easy. 

Up to now, we have seen that the main interaction tensor £,„ acts as the interaction term or the whole internal 
stress (Eq. (|22f) '). Now let us do the following statistics. The model introduced in Sec. ^Jcan be a crystal infinitely 
large in the three dimensions, or a big cut bulk of it. The bulk may be cut in different ways. For example, in Fig. 
2, the two cut bulks (only twenty-five cells each, due to the figure size limit) are slightly translated to each other 
and have some common MD particles with the rest MD particles being image particles to each other. Let us take an 
unweighted average of the dynamical equation, Eq. ||SJ|, over all possible cuts translated no more than one period to 
each other in all directions. In this averaging, only the term Fl^_r changes to 



•s, 



(34) 



as the same P is assumed acting on surfaces of all bulks. 
So the corresponding equations in Sec. IIIII become 

MrTrc = -8+ ~P^ S, 
D ■ T = (^e" + P ) • , 

H' = (^e^ + • V , 

wu = (^£^ + • ^cT , 
wii = (^£^ -pj ■ ^(T^ , 

while Eq. H23() still holds due to Eq. Q for equilibrium states. 



(35) 
(36) 

(37) 

(38) 

(39) 



COLLISION BETWEEN PARTICLES AND WALLS 



The individual collision between a surface particle and a wall occurs instantaneously and "randomly" . As we 
consider average surface effect only, let us do statistics on these collisions. The velocity of particle i in cell T at the 
surface of the bulk is -I- T. As the bulk is expanding or contracting with periodicity, the wall at the surface T moves 
at speed T. As usual, suppose particles are evenly distributed in the bulk or cells, the collision rate between particle 
i of cell T and an infinitesimal surface area vector ds of the bulk is 



j^ti ■ ds, if Ti ■ ds > 0, 
0, if r,-ds< 0, 



(40) 
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where the direction of ds is from the inside to the outside of the bulk. For ■ ds > 0, only the particle's velocity 
component in the ds direction is reversed by the collision, i.e., the particle acquires a total momentum of 

-2m,^—^ds (41) 
|ds| 

in each collision. The net external force on the R part (Fig. 1) due to collisions is then 

^'e,r= i ^^5I"^^l''^l'5(^0'^s = ^^m,|r/ g^ds, (42) 

where 



g (Oi) = (cos 0i + I cos 6'i I ) cos 6*^ , (43) 
Vi ■ ds , , , , 

|r,| |ds| 

Again, this force is dependent on the shape of the macroscopic bulk, which should not affect the properties of the 
inner part. So let us take the unweighted average of g [Oi) over any bulk shape. We also need to define two states as 
conjugate states with particle velocities {ri, r2, • ■ • ,r„} and {— i"!, — i"2j ■ ■ • 7 — i"n} respectively being the only difference 
between them. Further averaging g iOi) over conjugate states with equal weight leads to 

1 z"^'" r 1 

5=^y^ d^j^ gie)sinede^-. (45) 
Then replacing g (Oi) with 5 , Eq. (|42|l becomes 

Fk«= 3j7E"^^l^«l'S = l7^-S, (46) 

i 

where we got 

(47) 

i 

which is actually a scalar. 

The relation between the tensor P in the previous section and the constant external stress tensor T is then 

i 

but this relation can not be used to "remove" T in the dynamical simulation, as P is not known and T is given 
and fixed instead. Similar to Eq. ((JJ), we have 



' E,R 



T -ds = T -S. (49) 

R,sf 



When C is applied onto the surface of the bulk, the R part acquires additional momentum of the same amount per 
unit time. The dynamical equation of the R part becomes 

MrVrc - ^ E 1^*1' = Fl^r+ T -S. (50) 



In order to have the form of the dynamical equations closer to the ones in the PRMD, let us regard the term 



— mi\vi\^ I -S as a kind of internal force and move it from the left to the right hand side of Eq. H5()|l . yielding 



(51) 
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where the kinetic-energy term of the internal stress 



1 



By a similar deduction as in the previous section, we reach the dynamical equation 

D • T = + T ) • ^ct" , 

with the internal stress defined as 

TT = e + r . 
The similar simplified dynamical equations become 

H' = (^tT + T ) • , 

and 

WfL = (^tT + ^) • ^cT . 

For the case of constant external pressure, 

-P T , 



then Eq. H56|l tm-ns into 



For free systems, 



the equations become 



Wli = 



P ■ cr 



T =0, 

MrTrc = ^e" -S, 
D ■ T = ^e" • 

H' = ^e" • ^ct" 



(52) 

(53) 
(54) 
(55) 
(56) 

(57) 
(58) 

(59) 

(60) 
(61) 
(62) 
(63) 



as no wall and no t 



VI. TRANSPORT OF MOMENTUM 

In Eq. (|47() . the tensor C reflects the mechanical, regular, or real forces in collisions at surfaces, but we have not 
found what forces the kinetic-energy term of the internal stress r in Eqs. (|52|l and H54|l correspond to inside the 
bulk. Recently, Zhou proved that it corresponds to no regular forces (2^. Haile deducted this term by considering 
the transport of momentum, that is momentum transported byparticles passing through some plane(s) even without 
collision or any other interactions, in Appendix B.2 of his bookjj|. It seems that there should be an additional "force" 
in the pure transport of momentum. As particles may pass through numerous planes, we need make the additional 
"force" clear. 

Firstly, let us consider the simplest case, in which only one particle p with mass m is moving without being 
acted on by any regular force, so it is moving with a fixed non-zero velocity v. The corresponding dynamical equation 
is 



f = mv = 0, 



(64) 
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in which we see no force needed although the particle is passing through a great number of planes. Newton's First 
and Second Laws are always satisfied, but in the material system description (2^, in which the system always contains 
the same particles. Newton's Second Law is also usually written in the form 

where f , p, and t are force, momentum, and time, respectively. 

Let us view the same physical process in a spacial system description ^231. Imagine A and B are two neighboring 
empty space boxes separated by a plane PP'. In a short time interval At, particle p passes thorough the plane PP' 
from A into B. Let us define spacial system Sa as the system containing all particles in box A, and define spacial 
system Sb as the system containing all particles in box B. In these systems, the total number of particles may vary 
with time. The momentum increase is tov during At for system Sb- According to Newton's Second Law in the form 
of Eq. H65|l , there must be an additional "force" 

~* mv 

(66) 

acting on system 5*^ during this period. What provides such an additional "force" (no regular force at all)? When we 

do the same thing for system Sa, we find that another additional "force" f = — f 7^ is needed to satisfy Newton's 
Second Law for the same period At. These two additional "forces" are only needed when particle p passes from one 
system into the other. So we can regard this pair of additional "forces" as the interaction between the systems Sa 

and Sb, which means that Sa acts on Sb with the additional "force" f and Sb acts on Sa with the additional 

"force" — f . Apparently, this interaction obeys Newton's Third Law. In this case, the additional "force" changes 
an equilibrium dynamical process (Eq. (jHH)) into a non-equilibrium one (Eq. H66(l ^. but keeps Newton's Second Law 
satisfied. 

It is interesting to notice that if every aspect remains the same except that particle p moves from B into A with 
velocity — v, the additional "forces" between the two systems remain exactly the same. But this kind of additional 
"force" becomes more complicated when particle p is replaced by a bulk, with plane PP' cutting the bulk into two 
parts and moving at another velocity. It may be generalized to a physical picture in which particles are pushed 
to move by additional "forces" or particles activate additional "forces" when they move through a series of spacial 
systems even without any regular force acting on them. 

An application example for the additional "force" is an ideal gas contained in a cuboid container in a dynamical 
equilibrium state. Let us use plane PP' parallel to one end surface of the container to cut the gas into L and R parts 
with S as the cross section vector. With PP' fixed, both L and R parts are spacial systems, and particles may move 
from one into the other. It is easy to deduce the additional "forces" between the L and R parts as 

7= ± -S, (67) 

where r is defined in Eq. (|52|l and with L ot R part as one cell. As the regular external forces for collision of the 
gas and the container at the corresponding both end surfaces are 

f = T -S, (68) 

both L and R parts are in a dynamical equilibrium state. If we use material systems, S'j^ for example, in which those 
and only those particles residing in R at the time when the cut is performed always belong to 5^ (although they may 
pass thorough PP' later), the system S'j^ is not in an equilibrium state, but always satisfies Newton's Second Law, 
and no additional "force" is needed. In this example, the additional "force" changes a non-equilibrium dynamical 
process into an equilibrium dynamical process, with Newton's Second Law always satisfied. 

In the material system description, there is no transport of momentum into or from the system and no need to 
consider the additional "force" , and Newton's Second Law has its original form. If the spacial system description is 
employed, the net rate of the transport of momentum into the system, M, should be added to the acceleration side 

of the original equation of Newton's Second Law, then the corresponding additional "force" F= M must be added 
to the other side of the equation. As a result, the additional "forces" do not change the physics, but maintain the 
same physics when the point of view is changed. For continuous media, gas, or other systems with changing mass, 
the spacial system description might be a good choice, as "particles" are usually difficult to trace. 

The collision case in the previous section may be regarded as two superimposed dynamical processes. One is the 
low temperature dynamics. The other is the dynamics of an ideal gas in a container with changing volume and shape. 
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For the later, Eq. (|5()|l is in the material system description with — J^i "^j I part of the momentum 

rate acquired by the R part in the collisions, while Eq. I|51|l is in the spacial system description with r -S added to 
both sides as the net rate of the transport of momentum into the R part and the additional "force" acting on the R 
part by the L part, respectively. 

VII. CALCULATION EXAMPLE 

Using Eqs. Q and (|38|l with W — Mceii, some equilibrium structures of crystalline cobalt under constant external 
stress were calculated with an annealing technique adopted for low temperature, for which similar calculation for 
constant pressure was carried out in 1997j23j. The Lennard- Jones (L-J) 12-6 potential, v^'^{r) — 4£((cr/r)^^ — (tr/r)^) 
with e = 3.1225 x 10~^°J and a = 2.3059A was employed. Under normal conditions, cobalt has an ideal hep structure, 
with 2 atoms in the primitive cell and crystal parameters ao — 2.514A and cq = 4.105A. Table 1 presents the calculated 
structures for three cases, with a kept along the positive x-axis, b kept in the x — y plane with a positive y component, 
and c kept with a positive z component. In Case 1, the data are reproduced for the constant pressure of 1.0 atm. 
Then an additional pressure of 4.0 x 10^ atms is applied along the negative z-axis direction onto the top surface of 
Case 1, and the results are shown as Case 2. We see that the structure is still hep, but with |c| much decreased 
and |a| = |b| significantly increased due to the strong pressure along the z-axis direction. Further additional stress 

components of P 1,3= P 3.1— —4.0 x 10** atms are applied to Case 2, shown as Case 3. As a result, the structure is 
no longer exactly hep: |a| 7^ |b|, and c is perpendicular neither to a nor to b. However the structure still shows some 
features resembling hep. So we have changed the cobalt's structure by applying external forces or stress. 

VIII. SUMMARY 

By keeping Newton's Second Law for MD particles and applying it to macroscopic half systems, we arrived at 
a group of coupled dynamical equations ( Eqs. and H53|) ') for periodic systems under constant external stress. 
They clearly show that the periods change dynamically and that the external forces also influence the motion of MD 
particles. 

Our dynamical equations also show that the system is driven by the imbalance between the internal and external 
stress. The internal stress also takes the summation form of the interaction term and the kinetic- energy term. We 
further developed the interaction term into the main interaction term (Eq. ^14^ ) and the rest (Eq. H33|l'). The former, 
reflecting the main interactions between the two halves of the system, can be easily extended to cases of more accurate 
forces beyond pair potentials. The latter, taking a fixed form in all cases, reflects the statistics of different choices of 
the system. The controversial kinetic-energy term was obtained with the statistics of the collisions between particles 
and walls, and that of the additional "force" in the pure transport of momentum applied. 

From the dynamical equations, it is easy to see that even for low temperature equilibrium problems, the periods 
can not be assigned fixed values arbitrarily. Actually, the periods and MD particle positions must be configured such 

that both Eq. 101 and Eq. ^ are satisfied. As the external stress P is symmetric, Em should also be symmetric. 

For systems with estimated periods, Eq. (|23|l can be used to calculate the external stress P needed to keep them in 
equilibrium states. For free systems, the calculated external stress must be zero. Also from the dynamical equations, 
this approach works with both primitive cells and super-cells, so super-cells and the minimum-image criterion are not 
necessary. The previous example was done in a primitive cell as MD cell and the structures were calculated. 

With simulations for a series values of constant external stress, this approach can also be used in studying system 
properties changing with constant external stress. 
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Tabic 1. Calculated structures of a cobalt crystal for three cases of external constant stress P , with cell edge 
vectors a kept along the positive^ ;r-axis. b kept in the x — y planc^ with a positive y component, and c kept with a 
positive z coiupoueut. 





Case 1 


Case 2 




(atms) 


( 1.0 0.0 0.0 ^ 
0.0 1.0 0.0 
0.0 0.0 1.0 j 




^ 1.0 0.0 0.0 ^ 

0.0 1.0 0.0 
^ 0.0 0.0 4.0 X 10^ j 




( 1.0 0.0 4.0 X 10^ \ 
0.0 1.0 0.0 

y 4.U X iU U.U 4.U X iU / 


|a| (A) 
|b| (A) 

IcI fA) 


2.5140 

2.5140 
4.1047 


2.6259 

2.6259 
3.7073 


2.6130 

2.6294 
3.7126 


0(a,b) 

e{h,c) 


60.000° 
90.000° 
90.000° 


60.000° 
90.000° 

90.000° 


60.206° 
91.133° 
90.563° 


r2 - ri 


0.66667a 
-0.33333b 
0.50000c 


-0.33333a 
-0.33333b 
0.50000c 


-0.33097a 
-0.33807b 
0.50000c 
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CAPTIONS 



Table 1. Calculated structures of a cobalt crystal for three cases of external constant stress P , with cell edge vectors 
a kept along the positive a;-axis, b kept in the x — y plane with a positive y component, and c kept with a 
positive z component. 



Figure 1. A sketch for the bulk of a crystal cut into left (L) and right (i?) parts by the plane PP' with cross section 
area vector S parallel to the surface area vector uh of a cell, where h = a, b, or c is one cell edge vector. The 
half-linc-ccU bar is made of the MD cell, cells h, 2h, 3h, . . ., till the right surface of the bulk. Newton's 

Second Law is applied to the R part for the dynamical equations of the periods. 

Figure 2. A sketch for two cut bulks of complete cells from the same crystal, with the center cells as their MD cells 
respectively. They are slightly translated and have some common MD particles with the rest MD particles being 
image particles to each other. Each bulk has its own PP' plane. The complete interaction term of the internal 
stress is obtained by taking an unweighted average of the dynamical equations over these cuts. 
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FIG. 1: A sketch for the bulk of a crystal cut into left (L) and right (R) parts by the plane PP' with cross section area vector 
S parallel to the surface area vector cth of a cell, where h = a, b, or c is one cell edge vector. The half-line-cell bar Bh is made 
of the MD cell, colls h, 2h, 3h, . . ., till the right surface of the bulk. Newton's Second Law is applied to the R part for the 
dynamical equations of the periods. 
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FIG. 2: A sketch for two cut bulks of complete cells from the same crystal, with the center cells as their MD cells respectively. 
They are slightly translated and have some common MD particles with the rest MD particles being image particles to each other. 
Each bulk has its own PP' plane. The complete interaction term of the internal stress is obtained by taking an unweighted 
average of the dynamical equations over these cuts. 



